function r = geodist_wr ( lon, lat )
% GEODIST_WR:  distance between lat/lon coords in meters..
%
% uses the geodetic inverse formula (geodesy by %  A.R. Clarke) 
% Actually, this is just a wrapper around geodist.  It agrees
% with what the "distance" routine in the matlab mapping toolbox
% thinks to within maybe 1 percent or so.  Good enough for
% government work.
%
% USAGE:  r = geodist_wr ( lon, lat );
%
% PARAMETERS:
% Input:
%     lon, lat:  
%         in degrees 
% Output
%     r:
%         distance (range) between points in meters


n = length(lon);
for j = 1:n-1
	r(j,1) = geodist ( [lat(j) lon(j) lat(j+1) lon(j+1)] );
end

return


function [d,az1,az2]=geodist(input)
%  GEODIST uses the geodetic inverse formula (geodesy by
%  A.R. Clarke) 
%  input is a vector in the form (slat slon elat elon) such that
%            slat, slon - station latitude and longitude (degrees)
%            elat, elon - shot latitude and longitude (degrees)
%  output is:  d - distance in meters
%              az1 - azimuth station from shot (degrees)
%              az2 - azimuth shot from station (degrees)
%  adapted from a fortran subroutine used in the GLIMPCE expt.
%  
%  D. Hutchinson  1 September 1994






able = 6378.2064;
bake = 6356.5838;
easy = .006768658;

slat=input(1)*pi/180  ;                          %convert to radians
slon=input(2)*pi/180;
elat=input(3)*pi/180;
elon=input(4)*pi/180;

rat=(bake*bake)/(able*able);
cslat=cos(slat);
sslat=sin(slat);
celat=cos(elat);
selat=sin(elat);
cd=cos(slon-elon);
sd=sin(slon-elon);
ens=able/sqrt(1.-easy*sslat*sslat);
ene=able/sqrt(1.-easy*selat*selat);
ene=ene/ens;
a=-cslat*sd;
b= rat*celat*sslat -cslat*selat*cd +easy*celat*selat*ene;
az1= atan2(a,b);
a = celat*sd;
b = rat*cslat*selat -celat*sslat*cd +easy*cslat*sslat/ene;
az2 = atan2(a,b);
fkon = celat*cd*ene - cslat;
a = celat*sd*ene;
b = selat*ene - sslat;
b = b*rat;
fkon = fkon*fkon + a*a +b*b;
fkon = sqrt(fkon);
b = sqrt(easy/(1-easy));
a = b*sslat;
b = b*cslat*cos(az2);
beff = 1 + b*b;
fkor = fkon*beff;
beff = 1/beff;
bach = (a*a - b*b)*beff;
beff = a*b*beff;
c = zeros(1,4);
c(1) = -.1875*beff*bach - (1/3)*beff*beff*beff;
c(2) = .0046875 + .0375*bach + .25*beff*beff;
c(3) = -.125*beff;
c(4) = .04166666667;
d = c(1);
for k = 2,4;
   d = d*fkor + c(k);
end

d = d*fkor*fkor + 1;
d = (fkon*ens*d)*1000;
az1 = az1*180/pi;
az2 = az2*180/pi;



